Aufgabe ist es aus der Chemotion Repository (https://www.chemotion-repository.net/welcome) alle bisher hinterlegten Reaktionen zu klassifizieren, um anschließend sagen zu können wie sich die Klassen zwischen Edukten und Produkten verändert haben. Hierzu wird im Folgendem das verwendete R-Script beschrieben. Es wurden Daten von insgesamt 624 Reaktionen aus der Chemotion Repository am 30.04.2021 um 9:32 Uhr heruntergeladen.

Zunächst wird ein R-Skript mit den benötigten Helfer-Funktionen eingelesen, in der benötigte Funktionen ausgelagert werden, um die Vignette übersichtlicher zu gestalten.

source("../Scripts/viz-helper.R")

Die benötigten Pakete müssen geladen werden:

library(readr)
library(magrittr)
library(classyfireR)
library(plyr)
library(dplyr)
library(data.table)
library(networkD3)
library(ggplot2)
library(plotly)

Einlesen der CSV-Datei mit den Reaktionen und herausfiltern der nicht benötigten Zeilen:

data <- read.csv("../Data/reaction_export_30.04_0932.csv")
df <- data %>% dplyr::filter(!(type=="")) %>% 
  dplyr::filter(!(type=="3 solvent")) %>% dplyr::filter(!(type=="2 reactant"))

Im nächsten Schritt wird der Datensatz so formatiert, dass nur die Reaktionen in denen aus ein Edukt ein Produkt entsteht übrig bleiben. Dieser Datensatz wird dann anschließend in zwei Dataframes gesplittet, einen mit den Edukten und einen mit den Produkten.

more_start_mats <- df %>% slice(starting_mats(type))
one_start_mat <- df %>% anti_join(., more_start_mats)
one_start_mats <- one_start_mat[-c(51:53,80:82,325:327), ] 

start_mats <- one_start_mats %>%  
  dplyr::filter(!(type=="4 product"))

products <- one_start_mats %>% 
  dplyr::filter(!(type=="1 starting mat"))

Anzahl der Edukte und Produkte der zu analysierenden Reaktionen vor der Klassifizierung. Die Anzahl ist gleich, da nur Reaktionen mit einem Edukt und einem Produkt untersucht werden. Die Anzahl ist daher gleich der Anzahl der Reaktionen.

nrow(start_mats)
## [1] 434
nrow(products)
## [1] 434

Im nächsten Schritt werden die InChiKeys verwendet, um eine Klassifizierung der Substanzen aus den Reaktionen durchzuführen. Die Klassifizierung erfolgt mit dem Paket “classyfireR” (https://github.com/aberHRML/classyfireR). Die Funktion “get_classification” erzeugt ein JSON-Objekt für jede Substanz. Aus diesem Objekt können zum Beispiel Metadaten und die hierarchischen Klassifizierungsstufen ausgelesen werden. Für Substanzen die nicht klassifiziert werden konnten wird ein Objekt vom Typ “NULL” erzeugt.

Classification_List_start_mats <- sapply(start_mats$InChI, get_classification)
head(Classification_List_start_mats,1)
## $`GEYOCULIXLDCMW-UHFFFAOYSA-N`
## ── ClassyFire Object ───────────────────────────────────── classyfireR v0.3.6 ── 
## Object Size: 10.9 Kb 
##  
## Info: 
## ● InChIKey=GEYOCULIXLDCMW-UHFFFAOYSA-N
##   
## ● NC1=CC=CC=C1N
##   
## ● Classification Version: 2.1
##   
## kingdom : Organic compounds
## └─superclass : Benzenoids
##   └─class : Benzene and substituted derivatives
##     └─subclass : Aniline and substituted anilines
Classification_List_products <- sapply(products$InChI, get_classification)
head(Classification_List_products,1)
## $`BMIMNRPAEPIYDN-UHFFFAOYSA-N`
## ── ClassyFire Object ───────────────────────────────────── classyfireR v0.3.6 ── 
## Object Size: 14.3 Kb 
##  
## Info: 
## ● InChIKey=BMIMNRPAEPIYDN-UHFFFAOYSA-N
##   
## ● CC1=NC2=CC=CC=C2NC1=O
##   
## ● Classification Version: 2.1
##   
## kingdom : Organic compounds
## └─superclass : Organoheterocyclic compounds
##   └─class : Diazanaphthalenes
##     └─subclass : Benzodiazines
##       └─level 5 : Quinoxalines

Das erhaltene JSON Objekt aus ClassyfierR muss so formatiert werden, dass am Ende ein Datenframe mit der hierarischen Klassifizierung der Substanzen ausgestattet ist. Hierzu wurde eine Funktion classyfire2df_e() bzw. classyfire2df_p() geschrieben, die als Eingabe die Classification_List und als Ausgabe einen Datenframe für die Edukte bzw. Produkte hat.

class <- vector("character", length(Classification_List_start_mats[]))
df_start_mats <- classyfire2df_e()
class2 <- vector("character", length(Classification_List_products[]))
df_products <- classyfire2df_p()

Anzahl der Reaktionen, in denen die Substanzklassifikation erfolgreich war (Substanzklasse konnte bei Edukt und Produkt ermittelt werden):

sum(df_products$"1"!="")
## [1] 389

Die Daten der Ergebnisse der Klassifizierung werden im nächsten Schritt so vorbereitet, dass am Ende ein Dataframe mit den Klassen der Edukten und Produkten mit der dazugehörigen Anzahl entsteht.

dt1 <- setDF(data.table(df_start_mats[,2]))
dt2 <- setDF(data.table(df_products[,2]))
dt1$V2 <- dt2$V1

same_class <- dt1[which(dt1$V1 == dt1$V2), ]
same_classes <-ddply(same_class,.(V1,V2),summarize, size=length(V1) )

different_class <- dt1[which(dt1$V1 != dt1$V2), ]
different_classes <-ddply(different_class,.(V1,V2),summarize, size=length(V1) )
colnames(different_classes) <- c("source", "target", "value")
head(different_classes,9)
##         source                                    target value
## 1   Benzenoids Lignans, neolignans and related compounds     2
## 2   Benzenoids           Lipids and lipid-like molecules     3
## 3   Benzenoids                Organic nitrogen compounds     4
## 4   Benzenoids                  Organic oxygen compounds     4
## 5   Benzenoids                   Organohalogen compounds     3
## 6   Benzenoids              Organoheterocyclic compounds    72
## 7   Benzenoids          Phenylpropanoids and polyketides    13
## 8 Hydrocarbons                                Benzenoids     2
## 9 Hydrocarbons                    Organosulfur compounds     4

Anzahl der Reaktionen in der sich die Substanzklasse nicht verändert bzw. verändert hat:

sum(same_classes$size)
## [1] 217
sum(different_classes$value)
## [1] 172

In der folgenden Abbildung wird die Größe des Datensatzes einmal mittels Barplot besser dargestellt:

reactions <- data.frame (value  = c(624,nrow(products),sum(df_products$"1"!=""),sum(same_classes$size),sum(different_classes$value)),
                  reactions = c("all","one educt/product","one educt/product (classified)","one educt/product (classified) same_class","one educt/product (classified) different_class"))

o <- ggplot(reactions,aes(reactions,value))+
  geom_col(fill="darkblue")+
  aes(stringr::str_wrap(reactions, 15), value)+
  xlab("reactions")
ggplotly(o)

Der vorbereitete Datenframe kann nun genutzt werden, um einen Sankey-Plot erstellen zu können. Hierfür werden dem Datenframe noch ID’s zu den Klassen der Edukte und Produkte zugeordnet, auf die die folgende Funktion aus dem paket networkD3 (https://christophergandrud.github.io/networkD3/) zugreift.

different_classes$target <- paste0(different_classes$target, '.')

nodes <- data.frame(
  name=c(as.character(different_classes$source), as.character(different_classes$target)) %>% unique()
)

different_classes$IDsource <- match(different_classes$source, nodes$name)-1
different_classes$IDtarget <- match(different_classes$target, nodes$name)-1


ColourScal ='d3.scaleOrdinal() .range(["#FDE725FF","#B4DE2CFF","#6DCD59FF","#35B779FF","#1F9E89FF","#26828EFF","#31688EFF","#3E4A89FF","#482878FF","#440154FF"])'

p <- sankeyNetwork(Links = different_classes, Nodes = nodes,
                   Source = "IDsource", Target = "IDtarget",
                   Value = "value", NodeID = "name",sinksRight=FALSE,colourScale=ColourScal,nodeWidth=40,fontSize=13)

In dem erstellten Sankey-Plot befinden sich auf der linken Seite die Substanzklassen der Edukte und auf der rechten Seite die Substanzklassen der Produkte, in denen sich die Edukte nach der Reaktion umgewandelt haben.

p

Sessioninfo:

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] plotly_4.9.3      ggplot2_3.3.3     networkD3_0.4     data.table_1.14.0
## [5] dplyr_1.0.5       plyr_1.8.6        classyfireR_0.3.6 magrittr_2.0.1   
## [9] readr_1.4.0      
## 
## loaded via a namespace (and not attached):
##  [1] clisymbols_1.2.0  tidyselect_1.1.0  xfun_0.22         bslib_0.2.4      
##  [5] purrr_0.3.4       colorspace_2.0-0  vctrs_0.3.6       generics_0.1.0   
##  [9] htmltools_0.5.1.1 viridisLite_0.3.0 yaml_2.2.1        utf8_1.2.1       
## [13] rlang_0.4.10      jquerylib_0.1.4   pillar_1.5.1      glue_1.4.2       
## [17] withr_2.4.1       DBI_1.1.1         lifecycle_1.0.0   stringr_1.4.0    
## [21] munsell_0.5.0     gtable_0.3.0      htmlwidgets_1.5.3 evaluate_0.14    
## [25] labeling_0.4.2    knitr_1.31        crosstalk_1.1.1   curl_4.3         
## [29] fansi_0.4.2       Rcpp_1.0.6        scales_1.1.1      jsonlite_1.7.2   
## [33] hms_1.0.0         digest_0.6.27     stringi_1.5.3     grid_4.0.5       
## [37] cli_2.3.1         tools_4.0.5       sass_0.3.1        lazyeval_0.2.2   
## [41] tibble_3.1.0      crayon_1.4.1      tidyr_1.1.3       pkgconfig_2.0.3  
## [45] ellipsis_0.3.1    assertthat_0.2.1  rmarkdown_2.7     httr_1.4.2       
## [49] rstudioapi_0.13   R6_2.5.0          igraph_1.2.6      compiler_4.0.5